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ABSTRACT 

The aerodynamic performance of turbomachines is limited by flutter, choking, and stall, the last of which 
being a topic of this report. Viscous flow considerations are not only essential in determining loss coefficients, 
but also in ascertaining loadings and in providing a fundamental understanding of the flow. However, because 
of the complexity of flows in turbomachiues flow problems of model character only are accessible to numerical 
treatment. The focus of this report is on steady, incompressible, viscous, 2-d flows through a single infinite 
row of equidistant blades. In order to compute these cascade flows PrandtFs boundary layer concept enhanced 
by a strong interaction model is applied to the Navier-Stokes equations. Since this approach enables the 
integration of boundary layer equations in the presence of backflow, it is an attractive alternative to much 
more expensive Navier-Stokes solvers. 

Cebeci and collaborators, who developed such a viscous-inviscid interaction method for single airfoils, 
modified and extended their code to the cascade problem. The present contribution consists of a detailed 
description of the theory' and of conducting a few sample calculations of test-character. Flows past staggered 
cascades, using the ten percent thick airfoils of the NAC A-65-series compressor blade family as cross sections, 
were studied over a range of Reynolds numbers and air inlet angles. Difficulties were encountered in cases of 
low Reynolds number and highly loaded blades. The apparent success of the single airfoil code and its close 
agreement with experiments could not yet be extended to the cascade problem, partly because of the more 
demanding airfoil shapes used in turbomachiues and partly because of the inadequacy of the turbulence 
model at low Reynolds numbers. Improvements concerning both numerics and physics are necessary to turn 
the current cascade code into a robust and efficient tool to analyse cascade flows of practical interest. 
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1. INTRODUCTION AND GENERAL CONCEPTS 

Stall, the phenomenon at which lift decreases with an increase in angle of incidence, is one of the limiting 
factors in the efficient operation of turbomachines. The flow regime at stall is characterized by flow separation 
leading to regions of reversed flow, which do not only determine the loss behaviour but cause significant 
changes in loading and overall lift force. Flow separation is also present in another flow regime, Le. bubble 
transition, which is a low Reynolds number flow phenomenon, comprising laminar separation, transition, and 
turbulent reattachment. Because of economical requirements on turbomachines, which result in higher blade 
loadings, the aerodynamic community is forced to provide a better understanding of viscous effects in both 
attached and separated flow. The increasing costs of wind tunnel experiments and the decreasing computing 
costs encourage the development of new numerical toob. Solutions of the full Narier-S tokes equations are 
still reserved for supercomputers, and on the other hand, classical boundary layer codes are restricted to 
attached flows. Therefore, engineers have developed procedures, which do the job of a Navier-S tokes solver at 
a lower expense of computing time and storage. Such procedures are based on the vbcous-invbcid interaction 
concept which will be used in the present analysb. 

What are the ingredients of a successful viscous-in viscid interaction method? Following Prandtl's bound- 
ary layer concept the computational domain b divided into an invbcid region and a thin layer, close to the 
blade’s surface, where viscous effects cannot be neglected. The computational scheme solves separately for 
the flows in the viscous and in vise id domains, which was already a feature of classical boundary layer codes, 
but the essential point b the mechanbm how in rise id and viscous flows are matched. While classical bound- 
ary layer methods suffer from a pronounced hierarchy of riscous and in viscid regions, interaction models 
are supposed to exhibit the ability of mutual and immediate transfer of information between viscous and 
in vise id regions. In summary*, interaction methods must include an in viscid flow solution, a boundary layer 
method, able of handling both attached and separated flow, and some interaction modeL 

At this point we have to concede that our numerical possibilities are not even close to a solution of the 
flow inside a turbomachine. Unsteadiness, compressibility*, highly 3-d character, interaction of shocks and 
wakes, and complex flow boundaries, as found in turbomachiues. move a solution of the real problem far 
beyond today's potential. Because of the need to come up with some results, model problems are introduced, 
which retain only a limited number of the original features. Regarding geometry and decomposition of the 
3-d problem we adopt Wu's ideas to split the actual flow into two separate, but interacting 2-d flows, namely 
a blade-to-blade flow and a meridional-through flow. The blade-to-blade or cascade plane is defined as a 
surface of revolution formed by rotating a streamline of the meridional plane about the machine's axis. 
Such a projection of a single blade row. and our considerations are confined to single blade rows, reveals an 
infinite row of equidistant, similar, airfoil-shaped bodies, called cascade. Since we do not attempt a solution 
of the 3-d problem, the cascade problem will be isolated from the meridional-through flow. Further we 
do not allow any upstream disturbances, as generated by the wakes of upstream blade rows. The current 
calculations will also neglect the effect of viscous wakes emanating from the trailing edges of the considered 
blades. There Is no doubt that the flow in a turbomachine is highly unsteady and even the artificial cascade 
problem might experience unsteadiness due to a distorted upstream flowfield. However, because our model 
presumes steady flow, our ability to predict stalled flows will be somewhat restricted. Our limitation to 
incompressible flows is mainly related to the availability of efficient iuviscid cascade codes. This assumption 
dispenses us from taking shocks, shock reflections, and shock interactions into account. Finally the boundary 
layer approach imposes the restriction that viscous effects remain confined to a thin layer. Obviously, this 
assumption is violated in case of extensively separated flows. In view of all these drawbacks what advantages 
are being offered by iuviscid-viscous interaction methods and why is the interaction method superior to 
classical boundary layer methods and purely iuviscid flow solvers? The following levels of information can be 
associated with these procedures: Iuviscid flow solvers yield a pressure distribution and the overall lift force, 
which can be expected to be accurate only in case of attached and high Reynolds number flows. Boundary 
layer methods give additional information on skin friction distribution and overall drag force, but they fail if 
the flow separates. Interaction methods provide pressure and skin friction distributions, lift and drag forces 
for attached and mildly separated flows. 

The following sections will shed some light on the individual parts of this approach. The iuviscid flow 
equations are being solved by a standard panel method, which, mathematically speaking, provides a solution 
of Laplace's equation subject to Neumann type boundary conditions. Singularity methods, such as the panel 
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method, attempt a solution by an appropriate superposition of elementary solutions, named singularities, so 
that the boundary condition is satisfied. In two-dimensional flows only two types of singularities are required, 
namely sources (and sinks) to model displacement flows, and vortices to model lift. There remains the choice 
where and how to distribute these singularities. Smith and Hess [13] proposed to place them on the boundary 
surfaces which led to the first and still very successful panel codes. The first order approximation, which we 
are going to use, employs a representation of the airfoil by straight-line-elements, along which singularity 
strengths are assumed constant. The airfoil code of Smith and Hess was extended to cascades by Giesing [9], 
whose approach is virtually the same as the current one. except that Giesing’s method is derived and coded 
in complex variables, while the current one is done in real variables. Inviscid methods intended to interact 
with viscous flow solvers require some modifications. The commonly used flow tangency condition has to be 
replaced by the prescription of a wall transpiration velocity to make allowance for the displacement effect 
due to the existence of boundary layers. 

The second building block of the scheme is the viscous flow solution, being achieved by a finite dif- 
ference method. Finite difference methods are numerical tools to solve partial differential equations by 
discretizing the unknown functions at a limited number of gridpoints and subsequent rewriting of the differ- 
ential equations in form of algebraic equations, each of which applying to an elementary subdomain. The 
algorithm depends on the propagation of pertubations, such that elliptic equations require expansive field 
solutions, while parabolic equations permit marching procedures. Hence engineers try to avoid the elliptic 
Navier-Stokes equations and to focus on the parabolic boundary layer equations. An order of magnitude 
estimation, making use of the relation that the boundary layer thickness is much smaller than the streamwise 
length coordinate, justifies the reduction of the elliptic Navier-Stokes equations to the parabolic boundary 
layer equations. However, no gain without penalties: Only flow fields which obey the boundary layer ap- 
proximation ought to be treated, computations will break down at the point of zero skin friction due to the 
Goldstein singularity, and the downstream marching procedure will face numerical instabilities in regions of 
reversed flow. The first limitation Ls inherent to the approach. However, the Goldstein singularity is related 
to the boundary conditions. Catherall and Mangier [3] discovered that the boundary layer equations can be 
integrated through the point of separation, if displacement thickness or wall shear stress is prescribed in- 
stead of the external velocity. The third problem is tricky because of both up- and downstream propagating 
pertubations. which revert the problem locally to an elliptic one. The most common remedy, the so-called 
FLARE approximation [19|, drops the streamwise convective term in regions of backflow, thus getting rid 
of the numerical instabilities, but violating the physics whenever backflow velocities exceed small values. 
After having discussed the basic advantages and disadvantages of boundary layer methods, we continue with 
information on the currently employed method. Our viscous flow solution features Keller's box method [14|. 
which is implicit and second order accurate. In regions where separation is unlikely to occur the boundary 
layer equations are recast in terms of Falkner-Skan variables and solved subject to prescribed external veloci- 
ties. This direct method, which is applied in the immediate neighbourhood of the leading edge only, exhibits 
the classical hierarchical structure between viscous and inviscid flow regions, indicating that the outer flow 
imposes the pressure on the viscous layer. Soon after the leading edge or at the pressure peak, at the latest, 
the scheme switches from the direct to the interactive mode, which ensures the mutual interaction of inner 
and outer flow. While direct methods specify the pressure and inverse methods the displacement thickness, 
neither of these is prescribed in an interactive boundary condition. Both external velocity and displacement 
thickness are treated as unknowns and the boundary condition specifies a relation between the unknown 
streamfunction at the boundary layer edge, external velocity and displacement thickness. The concept of 
nearly constant boundary layer thickness, as provided by Falkner-Skan variables, has to be abandoned. The 
boundary layer equations are now written in terms of physical transformed coordinates and the Mechiil func- 
tion method by Cebeci accounts for the status of the external velocity as an unknown variable. Besides the 
crucial boundary condition at the outer edge of the viscous layer both methods employ the no-slip condition 
at the wall and the direct method requires initial conditions, which are generated by similarity solutions. 

In order to compute turbulent flows knowledge about Reynolds stresses, specifically turbulent shear 
stresses, is necessary. Unfortunately our understanding of turbulence is still quite inadequate. Therefore, 
we must use simple models based on hypotheses. The currently used model, the eddy viscosity concept by 
Cebeci and Smith [4], relates the turbulent shear stresses to mean flow properties on an empirical basis. This 
algebraic zero-equation closure allows to retain the same form of the differential equations and to apply the 
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same numerical procedures to both laminar and turbulent flows. In order to know where to start turbulent 
flow computations prediction of transition is necessary. Since again no exact methods are available, the best 
choice is to take experimental results. The attempt to predict transition is still in its infancy, resulting in 
rather limited applicability. Michel's method, which is used here, is based upon an empirical data correlation, 
narrowing considerably the range of applicability, while linear stability theories with empirical assumptions, 
like the e 9 -method offer some universality. Probably, transition and turbulence modelling represent the 
weakest members in the computational scheme. 1 " 

The interaction term describes the way how the inner viscous and outer inviscid regions communicate 
with each other. The transmitter of any information is obviously the edge boundary condition, which in case 
of a weak interaction works out as follows: The result of the inviscid flow computation, which is the pressure 
distribution, is imposed as boundary condition on the viscous flow computation. In return the result of the 
viscous flow solution, which is the displacement thickness, determines the effective shape of the body which 
is input to the inviscid flow solver. Basically, there are two possibilities to account for the displacement effect 
of a boundary layer on the outer inviscid flow: Either the blades are thickened by adding the displacement 
thickness to the blade coordinates or the flow tangency condition on the blade is substituted by a condition 
which prescribes a blowing velocity. Both way's can be incorporated either in a weak or strong interaction. A 
weak interaction will be sufficient if the viscous effect on the pressure is small permitting the above described, 
hierarchical manner of solving the complete flow problem. However, there are regions, like the trailing edge 
and separated flow, where viscous displacement might cause more severe changes in the local pressure field 
so that alternating treatment of external velocity and displacement thickness as in* and output appears to 
be too rigid a mechanism for the interaction of viscous and inviscid flow regimes. These sensitive areas call 
for equivalent and simultaneous treatment of displacement thickness and external velocity. Schemes, which 
meet these standards and hence do not possess a definite hierarchy, are named strong interaction models. 
Only they are capable of calculating boundary layers with flow separation. The approach taken here follows 
closely the ideas of Veldman [22|, who initially introduced the concept of strong interaction. 

The global organization of the scheme (see figure 1 on page 4) shows an iterative process in which 
alternately inviscid and viscous flow equations are being solved. Results are considered as converged if the 
external velocity distributions of viscous and inviscid flow solutions do coincide and if the results of the 
viscous flow solutions do remain steady over consecutive iteration cycles. Relaxation is a numerical tool 
either to ensure or to accelerate convergence. Viscous flow results are updated by relaxation just before they 
are transferred to the inviscid flow solver to serve there as boundary condition. The actually used relaxation 
formula was proposed by Kwon and Pletcher [15] and allows overrelaxation. 

The present report expands upon the work of Cebeci et aL [S], who developed a viscous-inviscid 
interaction method for incompressible cascade flows. The objective of the preseut work is to provide a 
complete description of the theory and to conduct a few sample calculations in order to begin the evaluation 
procedure. In chapter 2 a description of the panel method is given which is used to compute incompressible 
inviscid flows through arbitrarily staggered cascades of arbitrarily shaped blade sections. The solution of 
the boundary layer equations subject to an imposed pressure as well as to an interactive edge condition is 
discussed in chapter 3. which, additionally, reviews the used interaction and turbulence models. Flow results 
in form of pressure, skin friction, displacement thickness, and velocity profiles are presented and aualysed 
in chapter 4 for the NACA-65-series compressor blade family. The cascade configuration, namely a stagger 
angle of 45 degrees and a blade solidity of 0.78. was selected in consideration of available experimental data. 
Calculations were performed for Reynolds numbers of 245 thousand and 6 million and for air iulet angles 
ranging from 40 to 62 degrees. The report closes with a brief summary aud recommendations for future 
work. The reader is alerted to the fact that the preseut report is of an interim nature covering the period 
September 1985 through August 1986. 
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Fig. 1. Global organization of viscous-iuviscid interaction method 
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2. INVISCID FLOW METHOD 

Since the outer flow is assumed to be incompressible, a solution of Laplace’s equation will provide inviscid 
flow results. Laplace's equation will be subject to a Neumann type boundary condition, which prescribes 
a blowing velocity on the surface of the blades, resulting from the viscous displacement effect. In addition 
Kutta's condition, postulating zero load at the trailing edge, is imposed to establish a unique solution. 
The numerical implementation requires a couple of approximations which are associated with the term 
panel method. The following sections will cover a description of singularity methods in potential flows, the 
definition of the inviscid cascade problem, the concepts of influence coefficients and basic flows, and the 
derivation of farfield cascade relations enabling the statement of Kutta's condition for given inlet angle of 
attack and the calculation of inlet and exit flow conditions. 

2.1 Singularity methods In potential flows 

Irrotational flow fields [curlV = 3) permit that the velocity vector may be expressed as the gradient of a 
scalar function, called the velocity potential $, 



V = grad$. 

This definition together with the principle of conservation of mass for incompressible flows (div V* = 0) yields 
Laplace's equation, 

div[grad$) = 0 or A$ = 0 

which governs the steady flow field of incompressible inviscid fluids. The commonly used flow tangency 
condition is replaced here by a boundary condition that specifies a wall transpiration velocity 

— = V V f(*) on the surface of blades 

on 

Laplace's equation is a linear homogeneous second order partial differential equation, which is classified as 
elliptic, demanding a field solution. Two concepts have proven successful in the exact solution of Laplace's 
equation, the method of separation of variables and singularity methods. The latter ones achieve a solution 
by an appropriate superposition of elementary solutions, called singularities. This technique, underlying all 
panel codes, is based on Green's theorem and the principle of superposition, which applies to any linear 
partial differential equatiou 



If = 0. A = 0 A<&„ = 0 A$ ;V = 0 

then A$ = 0 with <5 = (7i$i + + . .. + C n $ n + . . . + C.v$.v 

where (7j, <7? C n , . .. , and Cy denote arbitrary' constants. Green's theorem, which is specific to Laplace's 

equation, states that the flow past arbitrarily shaped bodies may be represented by surface singularity 
distributions. This theorem simplifies the solution procedure considerably, because the flow at a field poiut 
can be computed by summing up the effects of all singularities instead of all field points. For 2-d flows the 
singularities to be used are sources and vortices, because the effect of displacement and thickness, respectively, 
can be handled by distributions of sources aud sinks, while the effect of lift is modelled by a distribution of 
vortices. The name singularity is attached to sources aud vortices, because they do satisfy Laplace's equation 
everywhere except at a singular point. 

Since the overall flow will be built up by simple flows, we review first the velocity potentials of basic 
flows which constitute the elements of our flow synthesis. The velocity potential of uniform parallel flow 
that is inclined with respect to the positive x-axis by au angle a is given by 

j?( jt. y) = ( x cos o + y sin a ) 

where U.^ deuotes the magnitude of the constant velocity vector. Radial streamlines aud concentric circular 
equipotentials characterize the flow field of a single source 

In \A 2 + y 2 

mis 
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where 9 denotes the source strength and the source is located at (0,0). Interchanging streamlines and 
equipotentials of a single source produces the flow field of a single vortex 

“y tl 

p(x, u) = - — arctan — 

2jt x 

where a positive 7 denotes the strength of a clockwise rotating vortex, located at (0,0). Our principal interest 
is in cascade flows so that our further considerations will deal with infinite rows of singularities. The velocity 
potential of an infinite row of equal and equidistant sources at the points(0,0), (0,±r), (0. ±2r), (0,±3r),... 
can be expressed in terms of hyperbolic and trigonometric functions 

E la\A 2 + (tf- kr)* = ^-lnf |( cosh cos 

*=+C O ^ * * 



A similar process leads to the velocity potential of an infinite row of equal and equidistant vortices at the 
points (0,0), (0,±r), (0, ±2r), (0, ±3r), ... 



<p(*,u) = - 



1_ 

2x 



fc=+oo 



E 

*=+00 



y- kr 
arctan 

■t 




2.2 Cascade flows 

Recalling that a cascade is defined as an infinite row of equidistant similar airfoil shaped bodies, we can 
conclude that the flow field past a cascade, which is inclined to a uniform steady onset flow, consists of 
identical subflowflelds. each of which is enclosed by two streamlines being displaced at the cascade’s spacing. 
Identical subflowflelds obviously can be caused only by identical singularity distributions on each body, so 
that the above defined infinite series will replace the expressions for single source and vortex, as used in the. 
airfoil formulae. On the basis of Laplace's equation's linearity we decompose the total velocity potential <5 
into that of uniform onset flow and the so-called disturbance potential j?. which accounts for the disturbance 
by the cascade of an originally uniform flow field 



$( x. y ) = £/*© ( x cos a + y sin a ) + x, y) 



where U ,? 0 and a are magnitude and direction (with respect to positive x-axis) of the vector mean velocity. 
(Note: The vector mean velocity is defined as the vector average of inlet and exit velocity vector). The 
disturbance velocity potential of a cascade . whose blades are arranged along the y-axis. can be written as 



*'•»>-/ IT 1 ” 



^cosh^-(/-/ 0 («)) - cos ^(y-yo(a))) 



/ tany(y-y 0 («)) 

- ^-arctan I * 



\tauli^*(x-Jo(^))> 



ds 



where the integration is carried out along the surface of the center blade, ( //o ) denotes the coordinates 
of the center blade, and r is the vertical spacing between adjacent blades. Since this integration cannot be 
performed analytically for arbitrarily shaped blades, the continues singularity distributions are discretized 
in the following fashion: 

1. Blades will be represented by straight line panels. 

2. The boundary condition will be satisfied at the midpoints of panels. 

3. Source strengths are assumed constant along each panel, but vary from panel to panel. Vortex strengths 
are constant over the blade. 

This unequal status of source and vortex strengths is due to the requirement that the number of unknowns 
must balance the number of equations. If the blade is approximated by N panels. .V boundary conditions 
applied at the midpoints of the N panels, plus Kutta’s condition at the trailing edge produce equations 
for N + 1 unknowns. Hence the N unknown source strenghts are determined such that the boundary 
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conditions are satisfied, while there is only one equation left to adjust one vortex strength. Introducing these 
approximations to the disturbance potential we get 



where 



rt +0/ J 

*?(*>!/) = S’ J ^ cosh ~{ x ~ x i~* cos 6 i) - cos “( V- Vi- 9 sin 0 j ) ) ] * 



tt 



7 " f ( tan-(y-yj-ssin^) \ 

- — XI / arctanj ^ ] <f* 

j=x y V tanh — ( x — x « cos OA ) 

—ij/2 ' r 



<r ; denotes the source strength of the j-th panel 
7 the vortex strength, which is constant for all panels 
{xj, jjj) are the coordinates of the midpoint of the j-th panel 

lj and Oj denote length of the j-th panel and the angle (positive counterclockwise) between surface 
direction of the j-th panel (positive clockwise) and positive x-axis. 



2.3 Concept of Influence coefficients 

Eventually we are concerned with velocities rather than the potential which is merely a vehicle to give 
information on the velocity field. Remembering that the velocity vector is the gradient of the potential we 
can apply this operator in any convenient reference coordinate system including local coordinate systems 
attached to panels. Now there are two velocity components we need to know: Firstly, normal velocities 
at panel-midpoints are required to satisfy the boundary conditions and secondly, tangential velocites are 
necessary to compute the pressure distribution on the blades. Both velocities will be computed by separately 
taking into account the influence of each panel resulting in 2 N contributions, one half induced by source 
distributions, the other half induced by vortex distributions. This procedure, called concept of influence 
coefficients, constitutes the numerical core of the panel method. Influence coefficients are hence defined 
as velocities induced by a unit singularity distribution of one panel Let »,• and U be outer normal and 
clockwise tangential directions at the midpoint (/,-, j/,) of the i-th paneL then influence coefficients due to 
source distributions take the form: 

Aiy represents the outer normal velocity at (/,-. y.) induced by a unit source distribution on the j-th panel 



A 0 = 







cos Oj) - cos — ( Ui— !jj— s sin Oj 




ds 



= - sin (Oi - Oj) V Xi . + cos (Oi - Oj) V yij . 



A tJ = 



1 

2 



for i j 
for i = j 



B t j represents the clockwise tangential velocity at ( x,*, Ui) due to a unit source distribution on the j-th panel 



+'j /’ 2 ( r 0 0 

S >J = £7 J ^7 1 ln [ cosh COS f)j ) - COS Ui - Uj - •» sin Oj) ) 

• * -',/2 ' ^ 

= cos(<9, - Oj) V XiJ - siu( Oi - Oj) V yiJ for / ^ j 
B, j = 0 for i = j 



ds 



with 






i 2 ^ ™ ) *2 ^ j) 

cosh - — cos — — 

cosh 2 irU, --^ +l) - cos 2 7 iklzLi±i \ 



Vyij = arctanAi + arctanA 2 - arctanA s | 
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■e the arguments Ai, Xj, A* are given by 

A*,) sing, + (fl-T,) cosA,]/, 

( AV) 1 + ( in - IV) 2 - Ml * - Ay) cos «y + 1 «- IV) »■» «,\ 

£i£izAi±ii cosh *K-.-V'i, sin *!«->>■) . r<»cl'atl sk |, ±i^Attd C03 ’Mzh 

_ r g g r j t 

• ffU.-AVn) 8 inh *(*,- a- /+ 1 ) £JJ g(tf,~r /+1 ) 1 cosh ^u.-A-y+x) :in g(w-rj 

3T X Tt It It 

*±izM cosb *'<-*/> sin djgrlil- -±«>zL). ^ cm *jjnzlA 

r(z,-A~yi 3inh «fci-.V/) , CM dflzjil + iiinlil C03t aiizAi siB ikizXA 

X XXX X X 



V, (J . and Vy.j denote x- and y-components (in a local coordinate system, whose x- and y-axes point in the 
clockwise tangential and outer normal directions of the j-th panel) of the velocity at (/,. y , ). induced by a 
unit source distribution of the j-th paneL Coordinates of boundary points (A/, I'/) and midpoints ( Xj-Uj) 
are related by 

*/ = j( Ay + Ay +1 ) = Ay + ^cos«y 

m “ jl>7 + ly+i) = iy + ^»“ 9 y 

Potential flow solutions frequently encounter troubles because of the multivalued arctangent. The following 
two provisions make sure that the arctangent, and subsequently the induced velocity, do not jump in value, 
where it is not supposed to do so: All arctangents must be evaluated in the range -i to +r and V yiJ . is 
calculated by addins up the effect of a single airfoil given by the first term, and the effect of a cascade 
without the center blade, represented by the second and third term. A minor difficulty, the computation 
of velocities at the midpoint of the inducing panel, can be overcome by Cauchy’s principal value technique, 
which permits integration of certain singular integrands. 

The influence coefficients due to a vortex distribution can be related to those due to a source distribution. 
Based on the fact that velocity fields of vortex and source are different by a clockwise rotation of 90 degrees 
only, the normal velocity due to a vortex equals the negative tangential velocity due to a source, and the 
tangential velocity due to a vortex equals the normal velocity due to a source 

< VORT EX _ nSOURCE 

A ij ~-Vij 

nVORTEX _ j SOURCE 

~ A iJ 



(Unless indicated otherwise influence coefficients will be those due to a source distribution.) 

2.4 Concept of basic flows 

Equipped with the concept of influence coefficients we can now deal with the satisfaction of the boundary 
condition, which is to determine source and vortex strengths such that velocities of given magnitude and 
direction be induced at panel midpoints: •• 



V * n Utf 1 T • \ 

„ = ( ► u’ 1 1 i 


for i = l... 


...V 


On, 






-V .V 






H A,]<r {0) - 7 J2 Bij— sin a - cos A.cos a ) 


for / = !,.. 


...V 



;=i ;=i 

To match the .Y + i unknowns (.V source and 1 vortex strength) by an equal number of equations the system is 
complemented by the empirical Kutta’s condition. The resulting system could be solved in a straightforward 
fashion by a linear equation algorithm: however, traditionally the flow is split into three basic flows, thus 
providing a very efficient solution for different angles of attack in case of no wall transpiration. Basic flows 
naturally divide into displacement and lift generating flows as follows: 
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1. Zero angle of attack flow (o = 0 deg, 7 = 0, wall transpiration): 

E + U*c sin B x - for « = 1, . . . , N 

2. Ninety degrees angle of attack flow (a = 90 deg, 7 = 0, no wall transpiration): 

N 

E A,.' 90 ’ = — t^cc cos Oj for i = 1, . . . , iV 

j=i 

3. Purely circulatory flow = 0, 7 = 1. no wall transpiration): 



tf TV 

E A^V 1 ’ = E Bij for i = 1 iV 

>=1 7=1 

Each of the above equations stands for a system of N equations with unknown source strengths only. Since 
all systems show up with the same coefficient matrix, the three sets of source strengths can be obtained by 
one execution of a linear equation solver with three different right hand sides. Although the concept of basic 
flows is primarily a numerical tool the following physical interpretation can be given: The first basic flow 
represents a pure displacement flow, whereby blades are inclined at zero angle of attack ( angle with respect 
to global x-axis. not blade chordline). The boundary condition of the first basic flow consists in a prescription 
of a wall transpiration velocity. The second basic flow is also a pure displacement flow, but at an angle of 
attack of ninety degrees and subject to a no-penetration condition. The third basic flow represents a purely 
circulator)' flow. Since a constant vortex distribution alone would not obey the flow tangency condition, 
normal velocities induced by this constant vortex distribution are balanced by an extra source distribution. 
In summary: Each of the basic flows satisfies some kind of boundary condition and each of them violates 
Kutta's condition. An appropriate linear combination of all three basic flows will yield flow at the required 
angle of attack subject to Kutta's condition: 

Total flow = “Zero angle of attack flow" * cos a 

+ "Ninety degrees angle of attack flows" * sin a 
+ “Purely circulator)* flow" * 7 

The vortex strength 7 is still unknown in the above equation, but provided it is known, total normal and 
tangential velocities at the midpoint of the i-th panel take the following form 






’ifUi-i/.) 



diij 

dflU, - !)i ) 
dti 



.V V 

U sin( «-<?<)+ E A ; jffj - 7 E Bi j 
j=i j = 1 

.v .v 

U^. cos( a - Bi)+ E Bij<r } + 7 E A,, 

j = 1 7=1 



with <Tj = <7j 0) cos or + <Tj 901 sin cv + <r. ^7 . The deviation of the outer normal velocities !*„ from the prescribed 
wall transpiration velocities \ t is a measure of the accuracy of the linear equation solver. 

2.5 Farfleld cascade relations 

For farfleld considerations the flow field past a cascade can be simulated by an infinite row of discrete vortices. 
Such an approach can be justified because the displacement-disturbances decay with increasing up- and 
downstream distance and because there remains a finite effect due to circulation in the farfleld. Using the 
formerly defined velocity potential of an infinite row of equal and equidistant vortices, this simplified model 
produces information about the horizontal and vertical disturbance velocities in the farfleld of a cascade 



= Hz 
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dx 



r 

— — arctau 
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where the small letters u and v denote the components of disturbance velocities, and the subscripts 1 and 
2 indicate streamwise positions at up* and downstream infinity. Magnitudes and angles of total inlet and 
exit velocities can be obtained by superposition of the above disturbance fiowfield and a uniform paralell 
flowfleld, which is inclined by an angle of a 



Abs( Inlet velocity): 



Air inlet angle: 



Abs( Outlet velocity): 
Air outlet angle: 



\Jv* + Vj 2 = \/(U, z ocosc»+ tii) 2 -t- ( 17^ sin or H- v x ) 3 



-fa+{h)'+ 



r \! sin o r 



tana, . 2. . t^a + r/ttf) 

Ui coso 



v't'J + v? - )Jtri + ( £)’ - 

<7 2 UccCosa 



where the capital letters U and V' designate the components of total velocities and U.^ is the magnitude of 
the vector mean ■velocity. (For an illustration of angles and velocity components see figure 2 on page 12). 

The calculation of forces requires either a pressure integration along the blade or the application of the 
momentum equation in integral form to an appropriate control surface. The latter method is to be preferred 
here because it involves fewer numerical steps and is said to be more accurate for small numbers of panels. 
Assuming steady 2-d flows of incompressible fluids the momentum equation can be written according to 



F 






pf v[vn) 

C 



ds 




pn ds 



where C denotes the control surface consisting of two vertical planes in the up- and downstream farfield and 
two similar streamlines, which are arranged such that they enclose a single blade. Making use of Bernoulli's 
equation magnitude and direction of the resultant gross force per unit span become 

f = \Jn77j = P u „ r with = \Ju? + ([ \\ + v 2 »/ 2 ) 2 

€ = arctanfF^/FJ =s a + 90 degrees 



The gross force turns out to act at a right angle relative to vector mean velocity, while no force is exerted 
in the direction of vector mean velocity. Lift is defined as the force perpendicular to vector mean velocity . 
Normalizing this force by a reference area and the dynamic pressure, based either on vector mean or inlet 
velocity, gives the dimensionless lift coefficient 






F _ 2U^T 

*ref(pU? ef /2) ~ CrefU? ef 



2.6 Kutta’s condition 

Kutta's condition assures that flows on the upper and lower surfaces merge smoothly at the trailing edge by 
postulating equal pressures and subsequent equal velocities there. Since the panel method does not permit 
the evaluation of velocities at boundary points. Kutta’s condition is satisfied in an approximate fashion at 
the control points of the first and last panel 



^ = -tf.v 



This condition determines the vortex strength 7, whose knowledge is necessary to synthesize the combined 
flow. Hence the scheme makes use of Kutta's condition after completing the basic flow solution and before 
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starting the combined flow computation. Cascade flows offer several input possibilities, two of which are 
being considered here. For a given average onset flow angle Kutta’s condition is a linear equation in the 
unknown vortex strength 7 



Ucc |cos(a- 0i) +cos(or- ^ jv)] + £ [(Bx ; • + Bjv/)(cosor^ 0) + sina<r^ 0, )| 

S N 

= — 1 [ E [Bij+ Byj)<r^ + 5Z (A-ij + 

l j=l J i - 1 J 



For given inlet air angle two unknowns, average onset flow angle or and vortex strength 7, remain to be 
determined, resulting in the need for another relation. Kutta's condition and the equation for the inlet air 
angle constitute a linear system in the two unknowns 7/ cos a and tana, which can be solved by any suitable 
procedure 



[l/oolsin^i + sin^v) + £ (B\j + S;v i )<rj- 90) | tan a + \ £ (Bi> + By j)<T^ + £ [Aij + -4.,v>)| — - — 

l jf= 1 J L ;=1 jf= 1 J cos a 

N 

= -U. x [ cosfli+cosMl- E (Biy + B[fj)<r[ 0) 

7=1 7 



r/ . ~ 

C^oc tan a + tan ax 

2r cosa 



One possibility to increase accuracy, especially if there are lower limits on the panel length at the trailing 
edge, is extrapolation of velocities at the trailing edge. This approach avoids the problem regarding singular 
behaviour of velocities at panel boundary points, and allows to impose Kutta’s condition exactly at the 
trailing edge. In case of linear extrapolation Kutta’s condition can be written as 



Vt % + 



V f| - Vi 



h + 



7 — - “ V '.v ” 

<2 



Iff + Av-i 



Av 



where the indices 1 and X refer to the first and last panel with numbering starting at the trailing edge on 
the lower surface and ending there on the upper surface. 

If the displacement effect due to the boundary layer is accounted for (nonzero wall transpiration velocity), 
it is advisable to evaluate the pressure distribution on the displacement surface itsell This step has to 
be taken because if the wall transpiration velocity exceeds small values, the assumption of no pressure 
gradient accross the ficticious inviscid flow region between blade and displacement surface ceases to be 
valid. The displacement surface is therefore the location, where total velocities are to be computed and 
where Kutta’s condition is to be met. While the scheme up to and including the calculation of basic flows 
remains unchanged, the scheme then proceeds with the computation of two more sets of influence coefficients. 
These so-called off-body influence coefficients represent normal and tangential velocities on the displacement 
surface induced by source and vortex distibutions of unit strength, residing on the original blade surface. 
The scheme makes use of these coefficients in determining circulation by imposing Kutta’s condition off- 
body at the trailing edge and computing total velocities. On the assumption of a rigorous (inviscid) theory 
the displacement surface will not be penetrated by fluid particles. However, because the computation of 
the blowing velocities is based upon a flat plate model, total velocities on the displacement surface have a 
tangential velocity component and a nonzero normal velocity component, which is smaller by an order of 
magnitude. The magnitude of the total velocities multiplied by the sign of the tangential component are 
transferred to the viscous flow solver. 
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Fig. 2. Cascade geometry and its vector diagrams 
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Fig. 3. Net rectangle for finite difference approximations 




sources are distributed 

Fig. 4. Visualization of the blowing velocity concept 
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8. VISCOUS FLOW METHOD 

In 1904 L. Prandtl gave birth to boundary layer theory, which closed the gap between the exact science of 
hydrodynamics and the more engineering-orientated hydraulics. Prandtl showed that the flow field past a 
body can be divided into two regions : 

1. a thin layer close to the surface of the body and the ensuing wake, where viscous effects have to be 
accounted for, and 

2. the remaining region, where inertia terms are dominant over viscous terms so that the latter may be 
neglected. 

These simplifications turn the elliptic Navier-Stokes equations into the parabolic boundary layer equations, 
which restrain the upstream propagation of disturbances. Numerically, this change of characteristics trans- 
forms the algorithm from a field procedure to a marching method, which integrates the boundary layer 
equations for given initial conditions step by step proceeding in the downstream direction. The differential 
equations are complemented by boundary conditions, whose choice significantly influences the applicability 
of each scheme: 

1. Direct boundary layer methods employ the so-called no-slip condition, requiring zero normal and zero 
tangential velocity at the surface, and a prescription of the external velocity at the edge of the boundary 
layer. Such procedures encounter the Goldstein singularity [10] at the point of separation, leading to a 
numerical breakdown. 

2. Boundary layer equations can be integrated through the point of separation, if the above edge boundary 
condition is replaced by a prescription of either displacement thickness or wall shear stress, while the 
no-slip condition is kept unchanged. However, such an inverse boundary layer method coupled with an 
inverse in viscid flow solution suffers from very poor convergence. 

3. Interactive boundary layer methods utilize an edge boundary condition which possesses the elliptic 
character of the outer flow. This so-called interaction law. which prescribes a linear combination of dis- 
placement thickness and external velocity, and the no-slip condition provide the required three boundary 
conditions for this scheme. 

The numerical approach features a finite difference method, which recasts the continuity and momentum 
equations as a system of linear algebraic equations. Keller's box method [14] approximates derivatives bv 
central differences and linearizes the nonlinear terms around the solution of the previous iteration. To 
stabilize the numerical algorithm in the presence of backflow the FLARE-approximation [19] is adopted in 
which the streamwise convection term is dropped in regions of backflow. The resulting set of difference 
equations is solved by Keller's block elimination method which takes advantage of the sparsely occupied 
matrix. 

In the following descriptions of both direct and interactive boundary layer methods, a derivation of 
the interaction law. and minimum information on the turbulence model will be given. The section "Direct 
Boundary Layer Method" will elucidate the more basic ingredients, like discretization of the flow field 
and Newton’s linearization, while the section on interactive boundary layer methods will emphasize its 
pecularities. 

3.1 Direct boundary layer method 

Direct methods are applicable in the absence of flow reversals, which limits their use in the current approach 
to the immediate neighbourhood of the leading edge. They are supposed to generate initial conditions at 
the stagnation point and to allow an efficient integration of the boundary layer equations around the leading 
edge. Generally speaking, direct methods appeal to the weak interaction problem in which viscous effects on 
the pressure are required to remain small. Although the paper addresses the strong interaction problem, all 
primary features of the finite difference method may as well be demonstrated for the direct method. To begin 
with we consider steady 2-d flows of incompressible fluids, described in a curvilinear coordinate system with 
/ directing along the blade’s surface and !j denoting the outer normal direction. The velocity components u 
and i’ shall be determined such that they satisfy anywhere in the flowfield the equations 

dn dr 



and at the boundaries of the flowfleld the conditions 



0: u{x, 0) = 0, r(.t, 0) = 0 

U=y e : u(x,y e ) = u e {x) 



where b denotes 1 + vt/u. The equations are referred to as boundary layer or thin shear layer equations, the 
conditions as no-slip and edge condition. In spite of the box method's requirement to reduce the governing 
equations to a system of first-order equations, we introduce for the present a streamfunction ( « = 34’ /dy 
and v = -dt/dx), which naturally decreases the number of dependent variables. Since streamfunctions by 
definition satisfy the continuity equation, only the momentum equation is left, which becomes a third order 
equation in t' T 



dy d 2 4’ 34? d 2 4’ due d 

dy dxdy dx dy 2 <<e dx + U dy dy 2 



Application of the Falkner-Skan transformation contributes to the removal of the singularity at x = 0 and 
maintains a nearly constant boundary layer thickness in terms of the new variable. In contrast to similarity 
concepts the dependency on two independent variables, x and q, is kept. The Falkner-Skan transformation 
scales the normal coordinate y and the streamfunction 4’ with reference to local external velocity and local 
streamwise coordinate 







v/“< 



vx 



4’(z,y) 



With prime denoting differentiation with respect to q the momentum equation and boundary conditions 
become 

i bry + ^±1 //" +m[i-(/'r 2 ] = *(/'!“- r jj) 

q = 0: f'[x< 0) =0. /(x, 0) = 0 



where //i = ( x/u t ) (du e /dx) is a dimensionless pressure gradient parameter. 

The box method involves four steps: reduction to a first order system, discretization in terms of central 
differences and two-point averages, linearization, and solution of the resulting linear system. The first step is 
accomplished by the introduction of two additional dependent variables. U and V', which converts the third 
order momentum equation into a system of three first order equations 



f = u 

U ' = V 



(M 



■V + 'JLt± fv + „ l{ i 



The system is complemented by the boundary conditions, which read in terms of the new variables 



n = 0: U[x. 0) = 0, /(A 0) = 0 

n = <U : = i 

The second step deals with the discretization of the continuous flow problem by considering instead of the 
continuos functions /. U . and V a set of discrete values of these flow properties. Let the solution domain 
0 < x < xj, 0 < // < if j be covered by a rectangular mesh 

xx = 0. x, = + ki with 2 < i < I 

q i = 0, q j - q j^i -h h j with 2 < j < J and qj = q t 
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Any flow property, whether it is a function or its derivative, will be expressed in terms of nodal values and 
coordinates of the network. Triplets of flow properties, consisting of streamfunction and its first and second 
derivative with respect to i), are assigned to each node of the network 

{/(*,,„,), cr ( V(x„ 9j )} = {/j, crj, vj) 

Because of the parabolic behaviour of the boundary layer equations the solution at a certain streamwfee 
position, say x,, depends solely on the solutions of upstream positions, say x,_i, x,_j, ... v while no down- 
stream influence has to be considered. The overall solution works in a step-by-step fashion, propagating in 
the downstream direction. The use of first order equations and the technique of central differences allow to 
reduce the domain of dependence from all upstream x-stations to the immediate preceding one. Hence one 
step of the solution procedure writes the governing equations for a column of net rectangles ( called boxes) 
residing in the subdomain 

x,'_i < x < x, and 0 < q < rjj 

and solves subsequently for the nodal values at the downstream face of the rectangular shaped subdomain. 
The x-station being currently solved holds therefore the superscript while the superscript "/—l" denotes 
the known flow properties of the adjacent upstream location. As indicated by the term central differences 
the equations are satisfied midway between the nodes: 

The two ordinary equations are centered about (x,, qy-i/a), and 
the partial equation is centered about ( x<_i/ 2 , JJy-i/j). 

The problem can now be rewritten in terms of finite differences (see figure 3 on page 12) 
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Boundary conditions undergo the same discretization process and take the form 



U{ = 0, f[ = 0 
= 1 

One might question whether this represents a solvable system in which the number of unknowns matches 
the number of equations. The subdomain under consideration consists of J - 1 net rectangles, the flow 
quantities of each being related by a momentum equation. Two auxiliary relations, which link dependent 
variables to their //-derivatives, refer to the downstream face of each net rectangle. Recalling the three 
boundary conditions we have a total of 3 J equations, which face J triplets of unknowns, each of the triplets 
corresponding to a node of the downstream side of the rectangular subdomain. Unfortunately the unknowns 
appear in nonlinear combinations, so that the solution involves an iterative procedure. The variables are 
linearized around their values of the preceding iteration (and in case of the first iteration around the solution 
of the adjacent upstream station) 
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where k denotes the iteration counter. If these Newton iterates are introduced and subsequently terms that 
are quadratic in 6Uj K , 6Vj’ K ) are dropped, a linear system in the unknowns fl/V*, 6U' } ’ K , 6Vj' K is 

obtained 
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where the coefficients abbreviate the following expressions 
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Together with the boundary conditions this system is repeatedly solved until the Newton iterates 6 /j ,K , 5Uj K , 

6Yj ,K are small enough to be neglected. As frequently experienced in discretized boundary value problems 
the resulting set of linearized finite difference equations displays a block-tridiagonal structure which permits 
a fast algorithm being applied to the linear equation solution. Block-tridiagonal matrices are composed of 
submatrices, named blocks, of which only those residing on main and both adjacent diagonals have nonzero 
entries. Equations and boundary conditions of the discretized boundary layer problem require a special 
arrangement to exhibit this block-tridiagonal structure. The system consists olJ sets, each of which include 
three equations: 

a. First set of equations (index j = 1) 

1. Wall boundary condition prescribing no penetration 

2. Wall boundary condition prescribing zero tangential velocity 

3. Second auxiliary relation (linking U' and V) of the first box 

b. Intermediate sets of equations (index 2 < j < J — 1) 

1. First auxiliary relation (linking /' and U ) of the ( J-l)-st box 

2. Momentum equation of the ( y — 1 ) -st box 

3. Second auxiliary relation (linking U* and V') of the /- th box 
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c. Last set of equations (index j = J) 

1. First auxiliary relation (linking /' and 17) of the last box 

2. Momentum equation of the last box 

3. Edge boundary condition. 



Boxes are numbered, the first referring to that at the wall, the last referring to that at the boundary 
layer edge. With {<5*'*} = {6f'- K ,8U'- K ,6V'j' K } T and {rj K } = {(ri)y K ,(ra)y K ,(r 5 )y*} T denoting vectors of 
the unknown Newton iterates and the known right hand sides, respectively, the linearized system can be 
expressed according to 
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where the blocks are 3 x 3-matrices being defined as 
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for 2 < ; < J - 1 
for 2 < j < J 
for 1 < j < J - 1 

1 -hjf 2 0 

{*sYj K i*s) r (3i ir 
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and where the right hand sides are obtained from 

mY = f;.;' - !‘r l 

J { momentum equation 

(r.»r = ur l - u j-~ l + 

(riJ'^-O. (r. 2 )i'" = 0. (fj)y K = 0 



for 2 < ; < J 
for 2 < j < J 
for 1 < j < J — 1 



The variety of indices explains itself by the 2-dimensionality and nonlinear character of the problem: 

1. The superscript indicates the location along the blade surface. 

2. The subscript ”/* denotes the location perpendicular to the blade surface. 

3. The superscript “k" designates the iteration counter. 

The unknowns are the Newton iterates of streamfuuctiou and its first and second derivative at the / -th 
streamwise position of the K-th iteration all across the boundary layer. The above system of linear equations 
can be solved most efficiently by Keller's block elimination method, which consists primarily of two types of 
recursions. A forward step eliminates the lower diagonal of submatrices, wheras the backward step solves 
the remaining system from bottom to top. 
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Note on discretization: The averaging of ail terms is carried through such that the number of generated 
terms is a minimum, which turns out to be trivial for single quantities, but has some impact on the product 
terms. Furthermore, averaging makes no difference whether acting on dependent variables or dependent 
parameters like "m" and n b". Let a and d be any dependent flow quantities, then discretization is performed 
on the basis of the following rules: 



1 . Evaluation of single terms: 

2 . Evaluation of differential quotients: 

3 . Evaluation of product terms: 

rather than: 



= j («}+ a< i-i + a T l + <*£) 
tda y-V* a' - a'~ l 
\dx) ~ ki 

/dav _ aj ~ ay-i 

Ul}/;-l/J hj 

(««*)/- 1/3= j[(«0y+(«0y-i] = + 

(ad)y_i/j = — ( Qjdj + ajdj—i + idj + Qj—i dj—i ) 



Although the second way of evaluating the product term is just as valid as the first one, the more straightfor- 
ward fashion of the first formula is preferred throughout the method. Care is advised in the reverse process, 
so 1/2 7*— ( acf) j— 1/2 because of the above agreement. 



3.2 Interactive boundary layer method 

The integration of boundary layer equations subject to an imposed pressure distribution terminates at the 
point of zero skin friction. Moreover, direct methods are restricted to weakly interacting regions, where 
viscous disturbances to the inviscid flow remain small Strong interaction arises for example from boundary 
layer separation or the rapid flow acceleration downstream of the trailing edge, both of which cause sub- 
stantial chauges in the external velocity distribution due to viscous displacement. Interactive methods are 
designed to permit integration through the point of zero skiu friction aud to account for strong interaction. 
The boundary layer approach Is maintained in strongly interacting regions, but accommodations apply to 
the status of external velocity and edge boundary condition. In contrast to direct and inverse methods both 
external velocity and displacement thickness are treated here as unknown quantities, reflecting the elliptic 
character of the outer flow. This introduces apparently one additional unknown into the viscous flow prob- 
lem. whose solution can be achieved either by an eigenvalue method or the Mechul function method, the 
latter of which being preferred here. The edge boundary condition of the direct problem is supplemented 
by the so-called interactive boundary condition, which relates the unknown external velocity with its iuvis- 
cid and '’displacement-perturbatiou-related" contributions. Boundary layer equations in the following form 
constitute a system in the unknown functions u(x. y), v(x , and u c ( x. y) 



dn dv 
r — H — = 0 
dx ay 

dn dn dn e d ..dn x 

11 T~+ " — = u e — +u—(b—) 
Ox dy dx dy dy 



0 = 



dn e 

dy 



The system consists of continuity equation, x-momentum equation, and the seemingly unnecessary y- 
momentum equation, wherein pressure in the momentum equations has been expressed in terms of the 
external velocity. The Mechul function approach assumes that the external velocity be a function of two 
arguments, resulting in the need for the trivial y-momentum equation. The reason for considering */*(/. y) 
rather than «*(/) is of a purely numerical nature. i.e., such a provision allows an easy setup of finite difference 
equations, avoiding nonlinear eigenvalue techniques. The governing equations are complemented by proper 
boundary conditions: The velocity components a and v are required to satisfy the no-slip condition at the 
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surface of the blade and the horizontal component must merge smoothly into the outer flow at the boundary 
layer edge 

y = 0: u{x, 0) = 0, v(z, 0) = 0 



U=U * : = u e (.r,0«), ««(•*, t/«) = «*/(*) + ~~ 



£ 



with u«/(^) denoting the inviscid velocity distribution and the second term, named Hilbert integral, approx* 
imating the perturbation velocity due to viscous effects. 

Interactive methods can be employed in both attached and separated flow regions .while direct methods 
cannot do so because of their breakdown related to the Goldstein singularity, nor inverse methods because of 
their poor convergence rates. Since interactive methods offer such a favourable behaviour, they are used in the 
current scheme for the main parts of the blade. Only the stagnation point singularity prohibits an exclusive 
use of the interactive method. The steps which turn the partial differential equations of the interactive 
problem into a linear system of algebraic equations resemble those of the direct method so that only major 
steps and their results will be repeated here. After the introduction of a streamfunction (u = dii'/dy aud 
t' = -dil> /dr), the equations undergo a transformation, which scales the normal coordinate y. streamfunction 
t, and the external velocity u e with reference to a constant velocity uo and the local streamwise coordinate 
x 



v) 




where u 0 is taken as the magnitude of vector mean velocity. The concept of constant boundary layer 
thickness, attained by Falkner-Skan variables (with « e as reference velocity), has to be abandoned because 
the external velocity is unknown in interactive calculations. However, provided that the integration does not 
start in the immediate neighbourhood of the stagnation point, growth of the boundary layer thickness can 
be kept limited. In terms of these so-called semi-transformed coordinates equations .written as a first order 
system by means of two additional dependent variables U and V', and boundary conditions take the form 



(bV)'+}/V + . rU 



f' = U 
U' = V' 
.dJV 
dr 
W' = 0 






n = 0: tf(j-.O) = 0, /(z.0) = 0 

n = n e : U(x,i} e ) = rc*(x. y e ) 






1 f d_ 

Mo !T J (l£ 




N'U. 0e)'le - /(£.'/«) 



dj 

r-Z 



The discretization of the flowfield follows closely the above outlined procedure of the direct method, covering 
the generation of an orthogonal grid and the introduction of central differences and two-point -averages. The 
overall solution proceeds in the downstream direction, accounting for downstream travelling disturbances 
only. In regions of backflow disturbances propagate against the direction of integration causing numerical 
instabilities. On the assumption that backflow velocities are comparatively small, a stable integration can 
be carried out by adopting the Flare approximation (Flugge-Lotz and Reyhner). in which the streamwise 
convection term udu/dr is set equal to zero in regions of backflow. With 



^>-1/2 



( 



1 

0 



if t'j-1/2 > 0 

if Uj_ < 0 
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designating an " on-off-switch” of the streamwise convection term, the finite difference equations of the 
interactive boundary layer problem become 



h. 



nz j Jh . = i {U . j+uU) 

EhEk.Lvt+v*,) 

hj 2 ' J } - 1 ’ 



= *.-1/2 



21 ^ 1 ) - (v-^ 

' V— 1/2 ' ^ j—l/2. 



- 1/2 /‘ - 



in - tt' 

— ^ ^- = 0 



Boundary conditions are expressed in terms of nodal values ,too, whereby the evaluation of the integral 
occuring in the interactive boundary condition involves an approximation in the fashion of the panel method 
approach, leading to 

U{ = 0 , f{ = 0 

U'j = W', W^gi + cuiWiu-f'j) 

where £,• and c,-,- denote a parameter and the diagonal element of the interaction matrix, resulting from 
a discrete approximation to the Hilbert integral. Averaging as well as centering is supposed to obey the 
principle that the number of generated terms approaches a minimum. This entails ordinary differential 
equations, like the y-momeutum equation, being centered about the middle of the downstream face and 
partial differential equations, like the x-momentum equation, being centered about the middle of the box. 

A balance of unknowns, which occur here as vectors in four components. {/*, U‘, Vj, ITj} 7 . confirms 
the principal solvability of the system: The 7 quadruplets of unknowns match with a total of 4 7 equations, 
including 2(7-1) auxiliary relations. (7-1) x-momentum and (7—1) y-momentum equations, each of which 
corresponding to oue of the (7 — 1) net rectangles, and 4 boundary conditions. After linearizing this system 
around the values of the preceding iteration (iteration counter "k- 1"), respectively around the solution of 
the adjacent upstream x-statiou in case of the first iteration, we arrive at a linear system in the Newton 
iterates Sfj K .SUj' K .^Vj- K . SW}-" 



/>/';* - 6 /& - f(w;- K + wfr) = /;_v l - + hjv£ 

h 
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supplemented by the two components of no-slip condition, edge and interactive boundary condition 
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Terms have been grouped such that known quantities reside on the right hand sides, while unknown quantities 
appear left of the equal sign. The abbreviated coefficients in the momentum equation are defined by 



“ hj + 2 k( W-V* •O-1/2) + 4 ' i 

1 -a 

(»«>/“ = - -fey + '-fir Us’ui - fi-U) + 1 /i-.' 1 

(»♦!}" - + V £W + 

/ \».<c 1 1 ■ tn-1 \ . it,'*** -1 

l $ *lj - 2 ki ' V i-*/2 + k ;-i/2/ + 4 V ;-1 

/ , \i.*e J '«— 1/2 rr>.M-l .1 

('*); jfc”^ ‘'/-i/2 

l a «lj ~ ~ U J-i °i-U » 

1 *, ' 

1. \i.K _*«-l/2 n-i.K—l 

0*), 



and the right hand side of the very equation reads 

Xj — i /2 f / n*^\ < **“ l _ / rr \ * ,>c 1 4. T /**** 1 f 4 * K “ l 4. V* - 1 /*<.«— 1 _ /*— 1 **«.«— 1 \ 

+ [l U J J-l/2 V C/ );-l/2 '^7— 1/2 + V ;-i/2 ' J-l/2 + V j-l/2 'J-l/2 'J- 1/2 V ;-l/2 j 



- { ltV| /‘- |tV|, /-‘» + l -uvy-.\ n - ^[(n.-n;:; /2 - wx-'x,- 



1/2 J j 



Since the overall procedure involves a repetitive linear pattern to approach the solution of the nonlinear 
system, the linear equation solver has to be as fast as possible. Fast algorithms take advantage of specific 
matrix strictures, one of which is the block-tridiagonal structure, which can be enforced in this method by 
the following arrangement of equations and boundary conditions: 

a. First set of equations (index j = 1) 

1. Wall boundary condition prescribing no penetration 

2. Wall boundary condition prescribing zero tangential velocity 

3. Second auxiliary relation (linking U 9 and V) of the first box 

4. Trivial y-niomentuni equation of the first box first box 

b. Intermediate sets of equations ( index 2 < j < J - 1) 

1. First auxiliary relation (linking /' and U ) of the (j-l)-st box 

2. Momentum equation of the ( J — 1 )-st box 

3. Second auxiliary relation (linking U ' and V') of the 

4. Trivial y-momentum equation of the /- th box j-tli box 

c. Last set of equations ( index j = J) 

1. First auxiliary relation (linking /' and U ) of the last box 

2. Momentum equation of the last box 

3. Interactive boundary condition 

4. Edge boundary condition. 
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Following these instructions equations and boundary conditions can be summarized in matrix-vector form, 
offering the desired structure 
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where { 6 1 :*} = and {rf } = {(n)f (r,)^} r denote four di- 

mensional vectors of the unknown Newton iterates and the known right hand sides, respectively. The blocks 
in the three diagonals of the above matrix have the dimension 4x4 and can be obtained from 
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The components of the right hand side vectors can be calculated from the following formulae 

for 2 < j < J 
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The numerical solution of the above system can a^ain be achieved by Keller s block elimination method, 
which works very much like a Gauss’s algorithm, but .firstly, matrices are eliminated instead of scalars, and. 
secondly, quite a few manipulations can be saved because of sparse occupation. 



3.3 Interaction model 

The coupling of the boundary layer and external inviscid flow is referred to as interaction model. Basically, 
the interaction consists of thickening the effective blade shape by viscous displacement, which in turn causes 
changes in the surface pressure. The exchange of information between viscous and inviscid regions takes 
place via the boundary conditions only, which is the specification of either an impermeable displacement 
surface or a nonzero wall transpiration at the original surface in case of inviscid flow, respectively the 
prescription of either pressure, or displacement thickness, or a linear combination of both in case of viscous 
flow. The interaction is called weak if viscous effects on pressure remain small while situations where 
viscous disturbances to the inviscid flow field are substantial demand the application of strong simultaneous 
interaction. Interaction models fall into four types, the first three of which provide loose coupling of viscous 
and inviscid regions: 

1. The direct interaction method combines a direct inviscid and a direct viscous flow solver. This classi- 
cal approach achieves a solution by iterative matching of boundary layer calculations with prescribed 
pressure and inviscid calculations with prescibed displacement thickness. The alternate treatment of 
pressure and displacement thickness leads to a local breakdown of the procedure when slight changes 
in pressure entail significant changes in displacement thickness. 

2. In inverse methods the roles of displacement thickness and pressure are interchanged. Hence the in- 
viscid flow equations determine the displacement thickness distribution, which is imposed as boundary 
condition on the viscous flow calculation. The results of which is the pessure, which closes the cycle 
by being input to the inviscid flow computation. The hierarchical manner of solving the complete flow 
problem excludes this model just as the previous one from handling strongly interacting regions. 

3. The semi-inverse interaction method is composed of a direct inviscid and an inverse viscous flow solver. 
Both parts of the scheme process the same input, Le. displacement thickness, and generate the same 
output, le. pressure. After each cycle an updated displacement thickness is relaxed based on the 
deviation of the two pressure distributions, which should coincide upon convergence. 

The existence of a definite hierarchy between boundary layer and outer flow in the above mentioned methods 
gives rise to a limited rate of feedback between both regions. Strong interaction models incorporate the outer 
flow somehow in the boundary layer calculations, for example by the following interaction law: 

4. Simultaneous (or strong) interaction methods solve the boundary layer equations subject to an interac- 
tion law as outer boundary condition, which retrieves the elliptic character of the outer flow. No definite 
assignment of displacement thickness and pressure can be made to viscous and inviscid region. Rather, 
both quantities are treated as unknowns, related by the interaction law. The procedure emphasizes 
simultaneous solution for both displacement thickness and pressure. 

In the following the currently used interaction law will be formulated in terms of displacement thickness 
aud external velocity, the latter being related to pressure by Bernoulli's equation. Both inviscid external 
velocity, resulting from a potential flow calculation around the original blade. « e /, and a correction term due 
to viscous displacement. « e j, are assumed to contribute to the total external velocity 



««( J*) = + « ts(*) 

The contribution due to viscous effects can be evaluated by means of the blowing velocity concept, which 
has its origin in a paper by Lighthill (1G|. He proved that the effect of boundary layers on the outer flow can 
be represented by a surface distribution of sources. Nonzero normal velocities at the surface of blades, being 
induced by this very source distribution, are intended to displace streamliues from the surface such that the 
virtual displacement surface becomes a streamline (see figure 4 on page 12) 

<#*(/) 

dx «*(/) 

where r(jr. <5*) denotes the vertical velocity component on the displacement surface. Since it is sufficient 
to approximate the correction term » e s. further deductions make use of the thin airfoil approach. Firstly, 
streamwise upper and lower surfaces of the blade will be considered flat plates in this context, implying 
that the blowing velocity equals half the local source streugth, and .secondly, the displacement thickness is 
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supposed to be so small that the horizontal velocity components of the inviscid flow do not vary across the 
boundary layer 

= V{X ' 0) = v{x ' ^ ~j 0 Ty dV= 

Unlike vertical components horizontal velocities are affected by all sources distributed along the supposedly 
flat surface. Eventually the viscous correction term is given by the well known thin airfoil integral, which is 
referred to here as Hilbert integral 








As indicated by the integration limits interaction is restricted to a finite region < * < x e on either 
streamwise upper or lower surface. The numerical implementation of the interaction law requires some 
discrete approximation of the above thin airfoil integral. Adopting the panel method approach, which 
concerns here a piecewise approximation of the continuous blowing velocity d(u e 6*)/dx to allow piecewise 
analytical integration, the integral can be written as a finite series 
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with [e,*| denoting a matrix of interaction coefficients. Recalling that the boundary layer calculation at a 
streamwise position involves iterations, the inviscid contribution can be included in the total external velocity 
of the previous Newton iteration (identified by leading to a generalized version of the interaction 

law 

K 

(«.)'■•* = K)'-*- 1 + E *k [(«,*•)*•* - (M*)^ -1 ] 

fcss l 

Since displacement thickness does not belong to the dependent variables, (tf*) 1 **, the displacement thickness 
of the streamwise position currently being solved, must be expressed in terms of dependent variables to make 
allowance for its unknown status. With (<5*) <# * being replaced by yj - n'j K /(«e) ,# * and after separating 
known and unknown terms, the interaction law takes its almost final form 






with 



9i = («e) i *' c_1 + E *,*[(«.*•)*•* - (M*)*'*" 1 ] - CuiU'Uj - 

A=1 



Implementing this relation necessitates a known right hand side, which can be evaluated only in an ap- 
proximate fashion, because the terms ( h^*)*’* are not known yet downstream of the current x-location. 
To ensure interaction also over this region, these terms are taken from the previous iteration updated by 
some relaxation formula. With cu- = dk y/vxjuo and g* = g*/u o designating dimensionless interaction 
coefficients and right hand side, respectively, the actually coded interaction law can be written in terms of 
semi-transformed coordinates 

tun*) + % 

This equation, being used as an outer boundary condition in the viscous flow computation, relates the 
unknown external velocity and the unknown displacement thickness of the i-th boundary layer station, 
whereby displacement thickness has been expressed by streamfunction and external velocity. Because of the 
elliptic character of the outer flow, which has been incorporated in the boundary layer via the interaction 
law. solution requires a global iteration, consisting of several viscous sweeps to be performed over both the 
streamwise upper and lower surface. 
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Note on the calculation of interaction coefficients: Let the continuous function "external velocity times 
displacement thickness”, denoted by D, be discretized at a finite number of streamwise positions. The thin 
airfoil integral then can be approximated by a finite series of weighted ”D's” at the very locations 

K 



I r£JL. z c ik D k 
*L X di n-i k % 



The weights are the interaction coefficients with the first subscript indicating the streamwise position, 
where the correction term u t f is to be evaluated, and the second indicating the location, whose effect is 
accounted for. With a properly interpolated D-function, integration can be carried out analytically piece 
by piece. Provided the point under consideration does not fall within the limits of integration. D will be 
interpolated linearly. Piecewise linear functions can be built up by overlapping triangular distributions, 
integration over which yields the coefficients whose k £ i, k £ i-1, k £ i+1 
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A linearly interpolated D-function would lead to singular integrals for the coefficients k = 1, k = i-1, and 
k a* i+1. Therefore D will be approximated by a polynomial of degree 2 in the interval JTi_i < Z < /i+i . 
Splitting again into overlapping distributions, which this time are parabolic, and applying Cauchy's principal 
value technique permits integration of singular integrands: 

The coefficient at the middle of the inducing source distribution is given by 
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As indicated above, the overall solution is approached in an iterative process, in which alternately 
viscous and inviscid flow equations are being solved (see also figure 1 on page 4 ): 

1. Calculate the external velocity distribution ««/ in an inviscid flow field by means of the panel method. 
To account for the blade-thickening due to viscous displacement, specify a nonzero normal velocity 
(blowing velocity) at the blade's surface. 

2 . March through the boundary layers of both streamwise upper and lower surface using ( over the main 
part of the blade) the interaction law as outer boundary condition. 

3 . Check convergence and quit if the criterion is satisfied. 

4 . If the convergence criterion is not met. prepare for another cycle. Update the product-term "external 
velocity times displacement thickness" on the basis of the deviation between inviscid and viscous external 
velocity distribution 

(M*r* +1 = (MT A [i + W ( - 1)] 

where A denotes the counter of global iterations. Further, compute the blowing velocity distribution, 
which serves as boundary condition for the inviscid flow solution 

(W 8 * 3J(“«0'- A+I 



and proceed with the first step. 

3.4 Turbulence modelling 

The presence of additional unknown shear stresses in turbulent flows requires modelling assumptions to 
balance the number of unknowns and equations. Eddy viscosity models, one of which is used in the present 
method, relate turbulent shear stresses to mean flow quantities on an empirical basis. They draw their 
versatility from the convenience of maintaining the same approach and numerical formulation for both 
laminar and turbulent flows. The currently employed model features a zero equation closure, which presumes 
that the length scale can be specified by an algebraic equation. The turbulent kinematic viscosity is defined 
by two separate formulae, one for the inner region being based on Van Driest’s approach ( 2 l|, and the ocher 
for the outer region being based on a velocity defect approach 
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Continuity of the turbulent kinematic viscosities is established by defining y e as the distance from the wall 
where expressions for inner and outer region do agree. Since transition is not an instantaneous process, 
an intermediate status of flow, called transitional flow, is assumed between the laminar and fully turbulent 
regions. This region is taken into account by introduction of an intermittency factor, which smears out the 
step-shaped change from kinematic to eddy viscosity 



7 tr = 1 - exp 



u* 



( ^tr U 



7 *r,r 



- 1.34 



(jp - 





where R (tr denotes the Reynolds number based on external velocity and streamwise location x fr at onset 
of transition. A value of 1200 was originally assigned to the empirical constant However, numerical 

experiments seem to indicate that low Reynolds number flows can be better modelled by values below 1200. 
The parameter a in the outer region formula is obtained from 
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where the noudimensional factor F denotes the ratio of total turbulence energy production to the shear* 
stress-related turbulence energy production, evaluated at the location of maximum shear stress 
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The ratio of the time-averaged quantities above can be approximated by a function of Rt = r„/(-«V) r 
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Since this turbulence model is not validated for separated flow, eddy viscosities in regions of backflow 
correspond to those of the adjacent upstream station at the same //-coordinate. 



Transition, if not available from other sources, currently is predicted by an empirical data correlation 
expressed in terms of the Reynolds numbers based on momentum thickness and streamwise coordinate at 
onset of transition 
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Should laminar transition occur upstream of the computed point of transition, then transition is redefined 
at the point of laminar separation. 



4. SAMPLE CALCULATIONS 

This chapter will focus on sample calculations reporting on the influence of turbulence model and loci of 
transition, as well as on the sensitivity of calculations to the choice of the relaxation parameter. Further 
calculations, referring to the variation of air inlet angle and Reynolds number, are intended to demonstrate 
the code’s abilities to predict the phenomena of low Reynolds number and high angle of attack flow regimes. 
No comparisons with experimental data will be presented in this report. However, the severe deviation of 
early calculations from experimental results, together with numerical breakdowns and convergence problems, 
revealed the need for a modification of the turbulence model in the low Reynolds number regime. Accord- 
ingly. the transition length has been reduced, resulting in fewer numerical troubles and improved agreement 
of computational and experimental results. Since the development of the code is not completed yet and 
further improvements are to be expected, the final evaluation, including the presentation of comparisons 
with experiments, is still to come. All the results describe flows past cascades, whose blades are staggered by 
an angle of 45 degrees and whose solidity amounts to 0.78. Cross sections of the NACA-65 series compressor 
blade family were selected in consideration of available experimental data. In the following, results of test 
case computations and a brief analysis of them will be given. 

Figures 5, 6. 7. and 8 show the convergence histories as a function of the relaxation parameter ^*. which 
has a major influence on the number of iterations required to obtain converged results. Relaxation is meant 
to ensure or to accelerate convergence. The lift coefficient based on inlet velocity is taken as convergence 
indicator, although other choices are equally possible. Figure 5 presents the convergence behaviour of a 
high Reynolds number flow at zero incidence. The geometrical input consists of 51 blade coordinates, taken 
from the original NACA publication. The code works with recomputed coordinates, which are cosine- 
distributed along the chord length. Subplot (a) refers to 51 of these recomputed coordinates, subplot (b) 
to 97. The upper plot of figure 5 displays a desirable convergence behaviour, wherein all curves approach 
the same lift coefficient after a finite number of iteration cycles. Maximum acceleration is achieved by a 
relaxation parameter of 2.0 or slightly higher. Smaller relaxation parameters delay convergence, higher 
relaxation parameters can lead to numerical breakdowns. The lower plot of figure 5 indicates a numerical 
instability, which seems to be related to the panel code. There is no evidence that these oscillations are due 
to unsteady loci of transition. The oscillations of the lift coefficients are stable in amplitude and the mean 
of all oscillations tends towards a single value. Since the amplitude of the oscillations amounts only to half 
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a percent of the lift coefficient, results might even be regarded as converged. Probably not the panel code 
itself, but the procedure how coordinates are recomputed is responsible for the instabilities. The currently 
used interpolation functions are - unlike spline functions • continuous in values of functions only, not in 
derivatives. If the number of recomputed coordinates strongly exceeds the number of input coordinates, 
then the used interpolation algorithm might give rise to irregular shapes in the leading and trailing edge 
areas, where points are packed because of the cosine-redistribution. 

The same cascade configuration at the same incidence, but at a much smaller Reynolds number leads to 
numerical divergence and numerical termination of the calculation, as can be seen from figure G. Subplots 
(a) and (b) refer again to 51 and 97 recomputed coordinates. Calculations have been performed up to 
40 iteration cycles, of which only the first 20 are shown. Some of the computations break down between 
cycle 20 and 40, the remaining computations keep changing, but do not get terminated within 40 cycles. 
Furthermore, computations for different relaxation parameters do not approach a common lift coefficient. 
Numerical breakdowns are either due to a failure of the Newton iterations at a particular streamwise position 
to converge, or due to negative external velocities as an artificial result of extensively separated boundary 
layers. These difficulties are related to unrealistically large regions of flow separation, which are caused by 
the turbulence model. The original CS-modeL being based on high Reynolds number data, overestimates 
the length of transitional flow for low Reynolds numbers, which in conjunction with bubble transition can 
trigger dramatic deviations from realistic flow results. 

Transitional flow in the CS-model is accounted for by an intermittency factor 7 < r 
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This factor, being a transient function with values 0 at onset of transition and 1 at fully turbulent flow, 
avoids an abrupt transition from laminar to turbulent flow by smoothing out the step-shaped change from 
kinematic to eddy viscosity. The length of transitional flow is mainly determined by the empirical constant 
G+f fr , which was set 1200 in the original CS-uiodel. This choice makes the transitional region of low Reynolds 
number flows spreading too far into the turbulent region. Consequently, low Reynolds number flows with 
a G lfr of 1200 are more likely to separate and less likely to reattach in the transitional region. Numerical 
experiments (see also figures 9 to 15) imply that low Reynolds number flows require smaller constants G^ fr . 
Figure 7 shows the results of the same low Reynolds number flows as in figure 6. using a slightly modified 
turbulence model, but still computing the point of transition. Subplot (a) corresponds to a G~ fr of 120. 
subplot (b) to a G lfr of 12. respectively. Calculations seem to converge now. however at retarded rates of 
convergence when compared with the high Reynolds number case. Most of the computations end up with 
oscillations of the lift coefficient, which originate here from moving transition points. Fixing transition on 
both upper and lower surface turns the oscillating lift histories of figure 7 into smooth curves, which are 
given in figure 8; Transition on the uppper surface is assumed at midchord, transition on the lower at 70 
percent of chord. Fixing transition supresses a degree of freedom, because of which only carefully selected 
locations might work out. The above chosen loci reflect the most upstream positions, where transition ever 
has been experienced within 20 cycles. Recapitulating, reasonable results for a low Reynolds number flow 
could be accomplished in the investigated case by the following provisions: 

1. The extension of transitional flow has to be narrowed by decreasing an empirical constant in the CS- 

model. 

2. The loci of transition have to be fixed at reasonable locations. 

The following 6 plots offer some insight into the convergence histories of the low Reynolds number flows, 
whose lift histories are given in figure C.a and 8.a. respectively. Figures 9. 11. and 13 refer to the unsuccessful 
calculation, using the original CS-model. computed transition, and a relaxation parameter of 2. Figures 10. 
12. and 14 display the results of an equivalent calculation, but now using the modified turbulence model with 
a G- lfr of 120. fixed transition at midchord on the upper and 70 percent of chord on the lower surface, and a 
relaxation parameter of 2. The histories of pressure coefficient, local skin friction coefficient, and displacement 
thickness are presented for the first 10 cycles. Figure 9 and 10 give the histories of pressure distributions, 
showing the developments from purely inviscid results (indicated by iteration number 0) to those altered 
increasingly by viscous displacement effects. Pressure distributions are rather insensitive compared to the 



more dramatic changes of the lift coefficients during the iteration process. With the exception of the trailing 
edge, pressure distributions keep looking the same, only flow separation leaves some marks - a sort of flat - 
on the pressure distributions. The longer flats in figure 9 suggest a more extended flow separation than in 
figure 10. 

Figures 11 and 12 , which show the developments of local skin friction distributions, prove this assertion. 
Two flow phenomena can be read from skin friction distributions: Flow separation occurs at the point of zero 
skin friction and transition is accompanied by a steep increase in skin friction. A closeup of the numerically 
unstable computation of figure 11 indicates that laminar flow separation takes place a few percent upstream 
of midchord on the upper surface, but turbulent reattachment, which closes the bubble usually not too 
far downstream of the point of separation, seems to be omitted. At the streamwise position, where fully 
turbulent flow is established, backflow dominates so that the boundary layer is not able to return to the wall. 
The lower surface exhibits laminar separation too, but there the flow does form a bubble. These results are 
considered unsatisfactory, because one would not expect a severe trailing edge separation, extending over 
half the chordlength, in cascades with only slightly cambered blade sections at zero incidence, even when 
exposed to low Reynolds numbers. The results of figure 12. which were obtained by using the modified 
turbulence model and fixed transition, fit much better both our expectations and experiments. Forced 
transition makes the results in some way less realistic, but forced transition is necessary to avoid oscillations. 
Transition is followed on both surfaces by a steep increase in skin friction, being caused by the fuller velocity 
profiles of turbulent boundary layers. The locations of transition, which had been chosen specifically to 
circumvent numerical instabilities, seem to suppress laminar separation and bubble transition. The flow 
remains completely attached on the lower surface, while it separates on the upper surface close to the 
trailing edge. This turbulent trailing edge separation can be noticed for the first time in cycle 3. and a few 
cycles afterwards a steady separation over 6 percent of chordlength has been established. 

Displacement thickness measures the deviation of viscous flow from in viscid flow. Le. it tells how far 
streamlines are shifted from the surface by viscous effects. Uneven displacement on upper and lower surface 
causes a lift loss compared to the inviscid flow case. Figure 13 and 14. which present the developments 
of displacement thickness distributions, depict the differences of displacement thickness on the upper and 
lower surface, which tend to uncamber the airfoil Displacement thickness at the trailing edge - like the 
lift coefficient - proves to be an excellent convergence indicator. The breakdown of calculations in figure 
13 is preceded by an impressive growth of displacement thickness at the trailing edge on the upper surface. 
Constant values since cycle 8 signal convergence of that case shown in figure 14. Both figures have marks, 
left by transition and separation. Transition is usually accompanied by slight drawback in the steady growth 
of displacement thickness, which can be seen clearly on both surfaces. Separation dramatically increases 
the growth of the displacement thickness, which is well demonstrated at the trailing edge on both upper 
surfaces. 

Figure 15 compares the results for two different airfoil sections of the NACA 65-series compressor blade 
family, both of which are arranged in a staggered cascade of relatively low solidity. The variations of lift 
coefficient based on inlet velocity and turning angle as a function of air inlet angle are given for both inviscid 
and converged viscous flow computations. The blade sections are the 10 percent thick NACA 65-410 and 
NACA 65-( 12)10. whose camber Ls indicated by the first number after the dash, denoting the design lift 
coefficient in tenths for the isolated airfoil. The high Reynolds number did allow calculations to be carried 
out with the original CS-turbulence model and having the points of transition computed. In spite of the 
high Reynolds number inviscid and viscous results are clearly different due to viscous displacement effects, 
which are more pronounced in case of the more cambered blade section. The beginning breakaway of the 
viscous curves from their inviscid counterparts, observed in the left subplot, suggests onset of stall at an air 
inlet angle of 62 degrees for the investigated cascade. Calculations for air inlet angles at and beyond stall 
could not be performed successfully. Therefore the plots were terminated before approaching angles of deep 
stall. 

The following three plots reveal the dependence of local flow quantities on air inlet angle for that cascade 
flow, whose lift variations are given in the left plot of figure 15. The variations of pressure, skin friction, and 
displacement thickness - all in dimensionless form - are shown in figure 16. 17 and 18. The distributions of 
pressure coefficients, based on the dynamic pressure of inlet velocity, provide information about the blade 
loadings. Pressure peaks result from highly accelerated flow around the leading edge, when the stagnation 
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point is relatively far away from the geometric leading edge. Only those pressure distributions, which 
correspond to air inlet angles around the design lift coefficient, can avoid velocity-peaks at the leading edge. 
The flattening of pressure distributions in the trailing edge area indicates separation on the upper surface 
for the highest incidence angles. None of the flows exhibits a laminar separation on either upper or lower 
surface, which anyway would not be expected for such a high Reynolds number. At air inlet angles of 60 and 
62 degrees turbulent trailing edge separation (without reattachment) does occur, extending over 4 and 20 
percent of chordlength, respectively. Transition from laminar to turbulent flow regime is accompanied by a 
step-shaped increase of skin friction, which permits the identification of transition in the skin friction plots. 
With increasing incidence transition is shifted on the upper surface towards the leading edge, while the lower 
surface experiences a reverse mow of transition. The development of displacement thickness distributions 
gives a good indication how much flow is disturbed by viscous effects. Generally, the perturbation increases 
with higher incidence. However, for positive incidence the lower surface shows no significant changes in the 
distributions of displacement thickness. It is the difference between the displacement thickness on upper and 
lower surface, that entails the loss of lift. Since the difference is more pronounced at high incidence, those 
flows depart more strongly from the inviscid ideal. 

Finally, the last four plots are intendend to promote the basic understanding of viscous flows past 
cascades. Velocity profiles along pressure and suction surface are shown for four representative flow regimes. 
The vertical scales of the velocity profiles had been enlarged to achieve a good resolution, allowing the 
illustration of marginal backflow. Figure 19 and 20 investigate the effect of the Reynolds number on viscous 
flow results. Flows at Reynolds numbers of 6 million and 245,000 were chosen to represent the high and low 
Reynolds number regimes at an air inlet angle, which corresponds roughly to the design lift coefficient. Both 
cases had been examined and depicted in terms of other flow quantities on previous occasions in this report: 
The lift histories of the high Reynolds number flow in figure 5.a, the lift histories and the developments of 
pressure, skin friction, and displacement thickness for the low Reynolds number flow in figures 8.a 10. 12. and 
14. The more obvious effects as the Reynolds number decreases are a more pronounced growth of boundary 
layers and a shift of transition points towards the trailing edge. The locations of transition in figure 20 
have been specified ( to avoid numerical instabilities); however, they do reflect the most upstream positions, 
where transition ever occured in the computation. The absence of a somehow expected transition bubble 
in figure 20 is due to that forced transition, which, if prescribed upstream of the laminar separation point, 
can eliminate the bubble entirely or reduce its extent. Furthermore, these figures demonstrate the different 
shapes of laminar, turbulent, and separating velocity profiles. The rapid flow acceleration immediately at 
the surface, as is typical for turbulent profiles, contrasts the gradual flow acceleration of laminar profiles. 

The influence of the air iulet augle on boundary layer results is investigated in figure 21 and 22. Air 
inlet angles of 48 and 58 degrees have been selected to show flow regimes of cascades at negative and positive 
incidence. As the air inlet angle increases, velocity profiles on upper and lower surface differ more aud 
more, resulting in strongly asymmetric displacement on the upper and lower surface. Transition on upper 
surface moves with increasing air inlet angle towards the leadiug edge, while on the lower surface the laminar 
flow region spreads with increasing air iulet angle, eventually leading to completely laminar flow, as can be 
seen from figure 22. Low Reynolds numbers aud small or negative incidence angles can give rise to bubble 
transition, which is illustrated in figure 21. Flow on upper surface remains laminar until midchord, separates 
then, becomes transitional, aud reattaches as a turbulent boundary layer. Strongly decelerating flows can 
even cause the separation of turbulent boundary layers. The increasing region of backflow on the upper 
surface at the trailing edge Is responsible for the loss of lift, as the augle of attack approaches stall. All four 
velocity profile plots give a good indication, how boundary layers grow, how they develop regions of reversed 
flow, aud also how the external velocity changes over the blade. 
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5. CONCLUSIONS 

This report presents a numerical method to analyse viscous incompressible flows past 2-d cascades. Cebeci 
and collaborators developed a viscous-inviscid interaction method, which accomplishes a complete flow* so- 
lution by treating alternately inviscid and viscous flow regimes. Inviscid flow equations are solved by means 
of a standard panel method, viscous flow equations are approximated .by boundary layer equations, whose 
solution is determined by a finite difference scheme. The essence of the method is an interactive boundary, 
condition, which forces viscous and inviscid flow to merge smoothly at the boundary layer edge. It is the 
numerical implementation of this condition, which, by prescribing a linear combination of extemaFvelocity 
and displacement thickness, avoids the difficulties incurred when either direct or inverse boundary layer 
methods are employed. The current method is — like Navier-S tokes-solvers — capable of computing viscous 
flows including flow separation, but in contrast to Navier-S tokes-solvers, at significantly lower computing 
costs. 

The potential aud capabilities of the method have been demonstrated for flows past cascades of moder- 
ately cambered blade sections over a range of Reynolds numbers and air inlet angles. The investigated flows 
involved both laminar separation bubbles and turbulent separation without reattachment. The procedure 
seems to work reliably for high Reynolds numbers and small incidence angle, but encounters difficulties in 
cases of low Reynolds numbers and highly loaded blades. Numerical experiments indicated the need for a 
modified turbulence model in low Reynolds number flows, in which the original model overestimated the 
length of transition. The location of transition evidently has a considerable influence on computed boundary 
layer results. Numerical instabilities had to be overcome partly by specifying transition, which might cause a 
deviation from the true flow* behaviour because of surpressed laminar separation bubbles. Since the method 
is not yet capable of handling highly cambered blade sections or incidence angles corresponding to stall and 
post-stall, further improvements are necessary to obtain a more universal code. Future work should attempt 
a more accurate modelling of both separating aud transitional flow regimes. In the near future progress may 
be made by 

1. Including the effects of viscous wakes, which should improve our ability to predict separated flow. 

2. Relaxing the FLARE-approximation by introducing an upwind differencing scheme for streamwise 
derivatives. 

3. Reconsidering some of the numerical features, like interpolation and smoothing algorithms. 

Ultimately, a satisfactory analysis of viscous cascade flows will require 

4. A better understanding of the physics of transition, eventually leading to a more reliable prediction of 
transition. 

5. A more advanced interaction model to resolve the difficulties at high incidence angles. 
ACKNOWLEDGMENTS 

This research was sponsored by the Naval Air Systems Command under Contract N-6227I-85-M-0437. The 
author wishes to express his appreciation to Prof. Platzer for his support and many helpful discussions. 
Prof Cebeci's help and the access to his codes are gratefully acknowledged. 



31 



Lift Coefficient Lift Coefficient 



0.576 



0.660 - 



CJ-2.6 




Cascade Data 

NACA 66-410 
Blade Solidity: 0.78 
Stagger Angle: 46 deg 
Air Inlet Angle: 54.3 deg 
Reynolds No.: 6 mill 



0.646 7 ( a ) 5 t Coordinates 



0.640 



-I L 



J 1 L 



J 1 L 



6 8 10 12 14 

Iteration Number 



16 



18 20 




Fig'. 5. Lift histories of a converging case. Flow past a staggered cascade at high Reynolds number 
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. G. Lift histories of a diverging - case. Flow past a staggered cascade at low Reynolds number 
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Fig. 7. Lift histories of a slightly oscillating case. Flow past a staggered cascade 
at low Reynolds number with a changed turbulence model 
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Fig. 3. Lift histories of a converging case. Flow past a staggered cascade at low 
■Reynolds number with a changed turbulence model and fixed transition 
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Fig. 9. Development of pressure distributions with number of iterations for a diverging case 
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Fig. 10. Development of pressure distributions with uumber of iterations for a converging case 






37 





Fig. 11. Development of local skin friction distribution with number of iterations for a converging case 



38 




Fig. 12. Development of local skin friction distribution with number of iterations for a converging case 
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Fig. 13. Development of displacement thickness distribution with number of iterations for a diverging case 
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Fig. 14. Development of displacement thickness distribution with number of iterations for a converging case 
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Fig. 15. Lift and turning angle as a function of air inlet angle for two different airfoil sections in cascade 



42 



30.0 1.2 r 1 30.0 









Fig. 16. Variation of pressure distribution with air inlet angle for a staggered cascade 
at high Reynolds number 
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Fig. 17. Variation of local skin friction distribution with air inlet augle for a stag- 
gered cascade at high Reynolds number 
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Fig. IS. Variation of displacement thickness distribution with air inlet angle for a 
staggered cascade at high Reynolds number 
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Scales Cascade Data: 




Symboli: TR... Transition 

SP... Separation 

Reynolds Number: 6000000 Air Inlet Angle: 54.3 deg 



Fig. 19. Velocity profiles on a N’ACA 65-410 airfoil in cascade at high Reynolds 
number and moderate angle of attack 
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Symbols: TR... Transition 

SP... Separation 

Reynolds Number: 245000 Air Inlet Angle: 54.3 deg 

Fig. 20. Velocity profiles on a MAC A 65-410 airfoil in cascade at low Reynolds num- 
ber and moderate angle of attack 
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Scales 



Cascade Data: 




Reynolds Number 500000 Air Inlet Angle: 48 deg 



Pig. 21. Velocity profiles ou a N AC A 65-410 airfoil in cascade at low Reynolds num- 
ber and low angle of attack 




SP... Separation 
RA... Raattachment 



Reynolds Number: 500000 Air Inlet Angle: 58 deg 

Fig. 22. Velocity profiles ou a NACA 65-410 airfoil iu cascade at low Reynolds uutn- 
ber aud high angle of attack 
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